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Abstract 

We present a numerical method for the computation of the con¬ 
formal map from unbounded multiply-connected domains onto lemnis¬ 
catic domains. For t-times connected domains the method requires 
solving £ boundary integral equations with the Neumann kernel. This 
can be done in 0(£ 2 n log n) operations, where n is the number of 
nodes in the discretization of each boundary component of the mul¬ 
tiply connected domain. As demonstrated by numerical examples, 
the method works for domains with close-to-touching boundaries, non- 
convex boundaries, piecewise smooth boundaries, and for domains of 
high connectivity. 
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1 Introduction 

In the theory of conformal mapping for multiply connected domains (open 
and connected sets) in the extended complex plane C = C U {oo}, there are 
several canonical domains onto which a given domain may be mapped. The 
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most commonly considered canonical domains are slit domains (see, e.g., 
Chapter VIII in the book of Nehari [39]), circular domains, and domains 
with polygonal boundary. Conformal maps onto these domains have been 
intensively studied, and several numerical methods for the computation of 
these maps have been proposed. Slit domains are considered, e.g., in EEl 
I91 I30H33] . circular domains in [3 [53 El], and Schwarz-Christoffel maps for 
polygonal domains in [ilf5[f7]l8l fT0HT2] . 

In this article we consider the numerical computation of the conformal 
map onto lemniscatic domains, which are another type of canonical domain. 
A lemniscatic domain is a domain of the form 

i 

C := jiu € C : JJ \ w — aj\ mj > r|, (1) 

3 = 1 

where a±,... , an 6 C are pairwise distinct, m\, ..., m,£, r > 0 are real num¬ 
bers, and the exponents satisfy 

i 

rn i = 1 • ( 2 ) 

3 = 1 

Lemniscatic domains where introduced by Walsh [ 41 ] , who proved that if /C is 
an Ctirnes connected domain with oo £ /C, there exists a lemniscatic domain 
C as in (JT]) and a conformal map <L> : AL —^ jC normalized by 4>(oo) = oo and 
^(oo) = 1; see Theorem [2J] below for the precise statement. The conformal 
map onto a lemniscatic domain is a direct generalization of the Riemann map 
for simply connected domains, for if t = 1 in (]T]) , then C is is the exterior of 
a disk. 

In addition to Walsh [13], the existence of the conformal map onto a 
lemniscatic domain was shown by Grunsky PME], Jenkins [20] and Lan¬ 
dau |26j. The last paper also contains an iteration method for computing 
<h, which, however, requires knowledge of the harmonic measure of parts of 
the boundary of the original domain /C. Recently, two of the present au¬ 
thors have investigated properties of this map and constructed some explicit 
examples in m- 

A remarkable feature of Walsh’s conformal map is that it allows a di¬ 
rect generalization of the classical Faber polynomials, which are defined for 
simply connected compact sets, to compact sets with several components. 
The resulting Faber-Walsh polynomials, introduced by Walsh in [45], are 
likely to prove useful, given the vast number of both theoretical and prac¬ 
tical applications of the classical Faber polynomials. For further details on 
Faber-Walsh polynomials we refer to the recent paper [ 143] . 
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Since the construction of conformal maps onto lemniscatic domains is in 
general nontrivial and only a few explicit examples are known, it is desirable 
to have a method for numerically computing such maps. In this paper 
we derive such a method and study it numerically. More precisely, given 
an £-times connected domain /C with a sufficiently smooth boundary our 
method computes the parameters defining the corresponding lemniscatic 
domain C as well as the boundary values of the conformal map <I> : AH —>• 
C. The method can be considered an extension of the approach described 
in [30H33] for the computation, in a unified way, of conformal maps onto all 
39 slit domains identified by Koebe in |23j. The method described in |30H33| 
requires solving a boundary integral equation with the generalized Neumann 
kernel. Using the Fast Multipole Method (FMM), the integral equation for 
multiply connected domains of connectivity l can be solved numerically in 
0(tn log n) operations where n is the number of nodes in the discretization 
of each boundary component [35 1[36j . The method presented in this article 
requires solving l boundary integral equations followed by solving a system 
of ln+l non-linear equations. The values of the conformal map for interior 
points can then be calculated using Cauchy’s integral formula. 

In recent years, several numerical methods have been proposed for com¬ 
puting the conformal map of multiply connected domains onto different 
types of canonical domains; see [Xlli Hl2l[28l[30l - [3il[36llT7l[i9] and the refer¬ 
ences cited therein. However, most of these numerical methods are limited 
to certain types of canonical domains or original domains. In comparison, 
the approach using the boundary integral equation with the generalized Neu¬ 
mann kernel can be used for a wide range of canonical domains. Moreover, 
it has been successfully applied to domains of very high connectivity, with 
piecewise smooth boundaries, with close-to-touching boundaries, and with 
complex geometry; see [3^1(361138] . 

This paper is organized as follows. In Section [2] we state Walsh’s exis¬ 
tence theorem for the conformal map onto lemniscatic domains. We then 
give the definition of the Neumann kernel in Section [3j In Section 0] we de¬ 
rive equations for the boundary values of the conformal map <f> and for the 
parameters of the lemniscatic domain C. In Section [5] we use these equations 
for the derivation of a numerical method for computing and £. Numerical 
examples with five different domains are presented in Section [6j Concluding 
remarks are given in Section [71 
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Figure 1: Illustration of Theorem 12.11 Original domain (left) and corre¬ 
sponding lemniscatic domain (right). The red dots are the centers ai, 02 , < 23 . 


2 The conformal map onto lemniscatic domains 


The following result is Walsh’s existence theorem from [44] (see also [ 26l 
Theorem 4]), which shows that lemniscatic domains are canonical domains 
for certain ^-times connected domains. 

Theorem 2.1. Let /C be an unbounded domain in C with 00 £ 1C, and let 
T = dK. consist of £ closed Jordan curves Ti,...,Tf. Then there exist a 
uniquely determined lemniscatic domain L of the form © and a uniquely 
determined bijective and conformal map 


4* : JC ->■ £ 


with 4>(z) 


z + O 



near infinity. 


(3) 


Further <f> extends to a continuous bijective map from /C = /C U T to C, and 
for each j = the image of Tj under contains the point aj in its 

interior. 

The number t > 0 is the transfinite diameter (i.e., the logarithmic ca¬ 
pacity) of the compact set C\/C. 

The uniqueness of the lemniscatic domain and the conformal map in 
Theorem O is forced by the normalization condition of near infinity 
expressed in ©. Note that if c £ C is any constant, then <f> c = + c maps 

K bijectively and conformally onto the translated lemniscatic domain C + c , 
and 4> c satisfies the normalization conditions <3? c (oo) = 00 and 4>(,(oo) = 1 . 

Theorem IQ is illustrated in Figure Q] for a domain 1C bounded by three 
Jordan curves; see Example 16.51 in Section [ 6 ] for details. 

We will need the following lemma. 
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Lemma 2.2. In the notation of Theorem \2.11 let 

t 

U(w) = Yl(w- aj ) m T 

3 =1 

so that CL = {w € C : \U(w)\ > r}. T/ien [/'(re) / 0 /or w £ dC = {w £ C : 

I^HI = r l- 

Proof. We show that the zeros of U' are the critical points of a Green’s 
function with pole at infinity, and that these are not located on dC. 

The function U is analytic but multi-valued in C\{ai,... ,ag}. Let 0 < 
p < t and consider the auxiliary function 

F{w) = log U(w) - log (p), 

which is analytic but not single-valued in C\{ai,..., an}. Then the harmonic 
functions 

G{w) = Re(F(u;)) = log \U(w)\ — log(p) and H(w) = Im (F(w)) 

are conjugates, i.e., G x = H y and G y = —H x . The derivative of F is given 
by 

-j-y = F' = — (F x — iF y ) = — (G x + iH x — i{G y + iH y )) = G x — iG y . 

Since G is real-valued, U'{wq) = 0 if and only if G x [w o) = 0 = G y {w o), i.e., 
the zeros of U' are the critical points of G. 

The function G is Green’s function with pole at infinity for the ^-tirnes 
connected domain {w £ C : |J7(u;)| > p}, see [45, P- 28]. Its contour lines 
(level curves) 

Ao- = {w G C : G(w) = log(fj)} = {w £ C : \U(w)\ = ap }, a > 1, 
have the following properties [ 451 p. 67]: 

1. For sufficiently small a > 1, A CT consists of I contours, each surrounding 
exactly one of the boundary contours (of {w £ C : |£/(ti;)| > p}). 

2 . A a grows with a, in the sense that if a < a', then A a is contained in 
the interior of A CT /. 

3. If a increases and A a crosses through an m-fold critical point of G , 
the number of components of A CT decreases by exactly m. 

Now, since C, is -Gtimes connected, dC = A r / p has I components. Therefore, 
no critical points of G can lie on dC (or interior to dC), which finishes the 
proof. □ 
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3 The Neumann kernel 


From now on we assume that K, is a given domain as in Theorem 12.II which 
additionally has a sufficiently smooth boundary T, oriented such that 1C is 
on the left of T. More precisely, we assume that each boundary curve Tj 
is parameterized by a 27r-periodic twice continuously differentiable complex 
function r)j(t) with non-vanishing first derivative 77 . 7 (f) = dr]j(t)/dt 7 ^ 0 for 
t 6 Jj = [0,27t], j = 1,2,... ,T (A dot always denotes the derivative with 
respect to the parameter t .) The total parameter domain J is the disjoint 
union of the t intervals J \,..., Jg, 


J ■— |_| 7? 

3 = 1 


|J{(t,j) : t G Jj}, 

3 = 1 


i.e., the elements of J are ordered pairs (t,j) where j is an auxiliary in¬ 
dex indicating which of the intervals contains the point t. We define a 
parametrization of the whole boundary T as the complex function 77 defined 
on J by 

V(t,j) ••= rij{t), teJj, j = (4) 

We shall assume for a given t that the auxiliary index j is known, so we 
replace the pair (f. j) in the left-hand side of ((4]) by t. Thus, the function 77 
in dU) is written as 

??i(f), t <E Ji, 
v(t) = < : 

m(t), t G j(. 

Let H be the space of all real-valued Holder continuous functions on the 
boundary T. In view of the smoothness of the parametrization 77 (f) of the 
boundary T, any function (j) € H can be interpreted via <f{f) := cj>(r](t)) as 
a 27r-periodic Holder continuous function of the parameter t on J, and vice 
versa. Henceforth, in this paper, we shall not distinguish between </>(i) and 

0(77(f)). 

We define the Neumann kernel N(s, t ) for (s, t) € J x J by 

v{t) 


N(s,t) := — Im 
7r 


77 (f) - 77 (s) 


It is a particular case of the generalized Neumann kernel considered in 
Similarly, the kernel M(s,t ) defined for (s,t) € J x J by 


M(s, f) := — Re 
7r 


77 (f) 


77 (f) - 77 ( 5 ) 
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is a particular case of the kernel M considered in [38]. The above kernels 
have been used in [311132] for computing the conformal map from unbounded 
multiply domains onto the canonical slit domains; see also |35H37j . The Neu¬ 
mann kernel also appears frequently in the integral equations for potential 
theory; see, e.g., |2l[TglTg], f25ll47] . 

Lemma 3.1 (see BSD- (a) The kernel N is continuous with 


N (t, t) = 7^ Im 


v(t) 

V (t)' 


(b) When s,t € Jj are in the same parameter interval Jj, then 

M(s,t ) = cot - + Mi(s,t) 

Zn Z 

with a continuous kernel M\ which takes on the diagonal the values 


Mi (t,t) = — Re 

Z7T 


v(t) 


We define integral operators N and M on the space H by 


N/i(s) ^ Jj N(s,t)n(t)dt, s G J, 

M/i(s) ^ L M(s,t)/j(t)dt, s € J. 

The integral operator N is a compact operator and the operator M is a 
singular operator. Both operators N and M are bounded on the space H 
and map H into itself [H]. Finally, any piecewise constant function v G H 
defined by 

v(t) = vj for t S Jj, 

with real constants v 3 for j = will be denoted by 


v(t) = (ui, t € J. 


4 The boundary values of $ and the parameters 
of jC 


In this section we derive equations for the boundary values of the conformal 
map and for the parameters of the lemniscatic domain C. These equations 
will be used for the numerical computation of and C in Section [5] below. 
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In the notation of Theorem 12.11 the function maps the boundary T of 
/C onto the boundary of the lemniscatic domain £, i.e., for t £ J, 


t 

ni*w«))-^r=r, 

j =i 

or, equivalently, 

e 

^2 rrij log\<f>(r](t)) - af = logr. (5) 

3 = 1 

To determine the parameters of the lemniscatic domain and the boundary 
values of d>, we fix for each j = 1,... ,£ an auxiliary point a.j in the interior 
of the curve Tj, and define the functions 

lj(t) := -l°g|r/(t) -cij\, teJ, j = (6) 


The first equation for the boundary values of will be obtained from 
the boundary values of the function 


t 

f ( z ) = ^rrijlog 
3 =1 


( ~ a 3 

V Z ~ a 3 


(7) 


where the branch of the logarithm with log(l) = 0 is chosen. The function 
/ is analytic in the domain /C with /(oo) = 0. By ([5]) , its boundary values 
are given by 

= l°gT + 7 (t) +in(t), teJ, (8) 

where 


l 

rf(t) ■=^2m j 'y j (t) 

3 = 1 

m(*) ; = Im f(v(t))- 


t 

^2 m 3 log \v{t) -OLj |, 

3 = 1 


(9) 


Although the functions y j are known from ((S|), the constant logr and the 
functions 7 and [i are not known a priori. We now show how the boundary 
values ([2D of / can be computed from the functions 7 j. A key ingredient is 
the following theorem from m- 

Theorem 4.1. For each function 7 j from ©, there exists a unique real¬ 
valued function [ij and a unique piecewise constant real-valued function hj = 
(h\j ,..., such that 

= lj(t) + hj(t) + i t € J, (10) 


8 



are boundary values of an analytic function fj in 1C with fj(oo) = 0. The 
function gj is the unique solution of the integral equation 


1^1 

1 

3 

■ft 

<0. 

II 

1 

s 

cr? 

(11) 

and the function hj is given by 


hj = [Mg,j — (I — N) 7 j]/2. 

(12) 


Theorem 14.11 allows to compute the functions hj and jij from the known 
function 7 j, so that the boundary values of fj are known. The following 
theorem shows that / = ^ J=1 nrijfj holds, by relating g and r to the known 
functions hj and fij. 

Theorem 4.2. Let f be the function from ED and let the notation be as in 
Theorem \f.l\ Then the function g and the constant logr in (HD are given 
by 

t 

h{t) = ^2m j n j (t), (13) 

j =i 

l 

log t = 'Y^m j h j {t), (14) 

3 = 1 

and we have f{z ) = J2j=i m jfj( z ) * n 
Proof. Define the auxiliary function g in 1C by 

t 

g(z) = f(z) m jfj{z ). 

3 = 1 

Then g is analytic in 1C with g{ oo) = 0, since each of the functions / and fj, 
j = 1,2,... has these properties. In view of © and (fTOl) the boundary 
values of g are given by 

tit 

g{r](t)) = log r + 7 (t) + i g{t) - ^ m j7? -(t) - ^ mjhj{t ) - i mjjij(t), 

3 =1 J=1 3 = 1 

which, by @, simplifies to 

g(g(t)) = log r — ^2 mj hj (t) + i( g(t) - ^mjHj{t) V (15) 
3=1 V 4=1 J 
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Thus 


l 

Re i9(v(t))] = logr nrijhjit ), 

3 = 1 

i.e., the real part of the boundary values of the analytic function g is a 
piecewise constant function. Since g{ oo) = 0, then g is the zero function [29l 
p. 165]. Hence, (ThU) and (fl3l) follow from (fl5l) . □ 

In order to compute the boundary values dHJ) of /, it remains to compute 
the numbers mi, m 2 ,..., mg. 

Theorem 4.3. Let the functions hj = (h±j ,..., hg,j ) be as in Theorem \f.l\ 
The unknown £ + 1 real constants mi,... ,m^,logr are the unique solution 
of the linear system 





'O' 


^ 1,1 

hi, 2 • 

• hi,e 

-T 


m 2 


0 


^ 2,1 

h2,2 ■ 

h2,e 

-l 

A 

mg 

= 

0 

, where A = 

he, 1 

CN^ 

... ^ 

• hg,e 

-l 


_l°g T_ 


1 


1 

1 ■ 

1 

0 _ 


(16) 


Proof. As shown above, the £ + 1 real constants mi,..., mg, r satisfy (USD- 
Since the functions hi ,..., hg are piecewise constant, it is easy to see that ([2]) 
and (|14j) can be written in the form of the linear algebraic system (1161) . 

We next show that A is nonsingular. Suppose that [<?i,... ,qg,c] T is a 
(real) solution of the homogeneous linear system 


Then we have 



l 

^2qjhj(t) = c, 
3 = 1 

t 

Y * = °- 

3= 1 


(17) 


(18a) 

(18b) 
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Let the auxiliary function f(z ) be defined by 

l 

/O) = 5^/iO) 

1= 1 

where the functions fj are as in Theorem 14. 11 Then / is analytic in 1C with 
/(oo) = 0, and its boundary values are given by 

ill 

qnj(t) + X! + 1 X] ^i (*)> 

l=i l=i l=i 

which by (|U|) and (I18al) can be written as 

e t 

= c - lo g I v{t) - atj\ + qjHj(t). 

3=1 3=1 

Define a function g(z ) on 1C by 


l 

£l(z) = f{z)+ ^2qj\og(z - aj). (19) 

3=1 

The function g(z) is analytic in 1C but is not necessarily single-valued. For 
large |z|, we have 

log (2 - ocj) = log [z (l - = log z T log (l - , 

which implies in view of (|18bl) that 


l It 

^2 i q j \og{z-a j ) = log(z)^2qj+Y2 q 3 l °S ( X 

j=i 1=1 1=1 


a, 


og^ 1 

l=i 


a, 


Hence, the second term in the right-hand side of (1191) vanishes at oo. Since 
/(oo) = 0, we have g{ oo) = 0. Let the real function u be defined for z £ /CUT 
by 

u(z) = Re 0(z). 


Then u is harmonic in 1C, and satisfies the Dirichlet boundary condition 


u = c on T, 


( 20 ) 
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i.e., u is constant on the boundary T. The Dirichlet problem (1201) has the 
unique solution u(z) = c for all z £ /C U T, so that the real part of g is 
constant. Then g is constant in 1C by the Cauchy-Riemann equations, and 
g( oo) = 0 shows that g(z) = 0 for all z G 1C, and, in particular, c = 0. 
Then, (1191) implies that 


i 

f(z) = — qj log {z — a.j) for all z € 1C, 

3 = 1 

which is impossible unless <?i = q 2 = '"=% = 0 (since the function on 
the left-hand side is single-valued and the function on the right-hand side 
is multi-valued). Thus, the homogeneous linear system (|17j) has only the 
trivial solution [gi,... ,qi,c\ 7 = 0, and A is nonsingular. □ 

By obtaining the real constants mi,... ,rri£, we can compute the func¬ 
tions 7 and g from @ and (fl3l) . By (| 8 j), the boundary values of the analytic 
function / from (j 2 j) are given by 

lo S ( ^|) - o. ° J ) =1 °g r + TW + ( 21 ) 

This is the first equation needed to compute the boundary values of the 
conformal map $. Note that in (1211) only $>(r)(t)) and the complex numbers 
a,j are unknown. We will need one more set of equations to determine these 
quantities. 

Lemma 4.4. The boundary values of the function <J> and the i constants 
a \,..., ai satisfy the i equations 

s/ J '° l! (llp?) ,|,|< + “ r “ i = 1 ’ J = ( 22 ) 

Proof. Since the functions 

ifjiz) = log j = 1,2,... ,£, (23) 

V z a j J 

where the branch of the logarithm with log(l) = 0 is chosen, are analytic in 
the domain 1C including the point at oo with ipj( oo) = 0, we have 

--r f ipj(ri)dr / = Res if>j(z) = - Res ( - ) , j = 1,2,..., t, (24) 

27T1 J p z=oo J 2=0 \Z J 
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see [2TJ pp. 107-108]. Let 1C be the image of the domain 1C under the mapping 
r and the functions (jj (z) be defined on 1C by 


9j( z ) = ifrj 


, j = 1 , 2 ,...,* 


Then gj is analytic in 1C with gj( 0) = 0 for j = 1, 2Hence 

Res -4 ipj ( - J = lim (-) = lim = g 3 (0), j = 1,2,...,> 
z= 0 Z z \Z J zs- 0 Z \Z J z-¥0 Z J 

To compute g'(0), we have from 

<5>{l/z) - aj 


(25) 


9j{z) = log 


= log 


z<b{l/z) — ajZ 

1 — ajZ 


1 / z — a j 

For small \z\, the normalization ([3]) of implies z$(l/z) = 1 + 0(z 2 ), and 

1 — CljZ + 0(z 2 ) 


9j( z ) = log 


1 — ajZ 


so that 


which implies 


-aj + 0(z) 


— Otn 


^ ^ 1 — CljZ + 0{z 2 ) 1 — OtjZ ’ 


9j (°) = °lj — dj , j = 1,2,. .. (26) 

The assertion of the lemma follows from ( 1751 ) . (|24|) . ( 1751 ) . and ( 1751 ) . □ 


5 The numerical computation of the conformal 
map $ and lemniscatic domain C 

In this section we discuss the numerical aspects of the computation of the 
conformal map $ and the lemniscatic domain C based on the results from 
Section [H 

We first consider the numerical solution of the boundary integral equa¬ 
tions from Theorem 14.11 and the computation of the parameters m\,, rri( 
and logr from Theorem 14.31 We then consider the computation of the 
boundary values of # and of a\, ..., ai by solving the equations ([TIT) and ( 1771 ) . 
Finally, we discuss the computation of the values of at interior points of 
1C by Cauchy’s integral formula. 
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5.1 Computation of the parameters mi,... ,mg, logr 

The £ boundary integral equations CEB can be solved accurately by the 
Nystrom method with the trapezoidal rule [21125] (see [301131 1 1 3 5 1 136] for more 
details). Let n be a given even positive integer. For k = 1,2, ...,£, each 
interval J& is discretized by the n equidistant nodes 

2 vr 

Sk, P = {p- 1 )— £4 P = 1,2,... ,n. 
n 

Hence, the total number of nodes in the parameter domain J is £n. We 
denote these nodes by ti, i = 1 , 2 ,.. ., £n, i.e., 

t(k—l)n+p ^ Ji h 1)2, ...,£, p 1,2,... ,n. (27) 

For domains with piecewise smooth boundaries, singularity subtraction m 
and the trapezoidal rule with a graded mesh [24] are used (see [25]). By 
discretizing the integral equation CEB by the Nystrom method with the 
trapezoidal rule, we obtain the £n x £n linear algebraic system (I — H)x = y; 
see, e.g., [36] equation (59)] for an explicit formula for this system. 

Since the integral operator N is compact, the only possible accumulation 
point of its eigenvalues is 0; see, e.g., [25l p. 40]. Moreover, as shown in [36, 
Corollary 2], the eigenvalues of I — N are contained in the interval (0,2] 
(see also [25, Theorem 10.21]). Consequently, for sufficiently large n, the 
eigenvalues of the discretized operator I — B are contained in (0, 2] and they 
cluster around 1. Numerical illustrations of this eigenvalue distribution are 
shown in m Figures 4-5]. 

We solve the discretized system (I — B)x = y using the (full) GMRES 
method SB- Each GMRES iteration requires one matrix-vector product 
with I — B. This product can be efficiently computed using the Fast Mul¬ 
tipole Method in just 0 (£n) operations [HICLE]- It was already observed 
in [36], that the number of GMRES iterations for obtaining a very good 
approximation of the exact solution (relative residual norm < 10~ 12 ) is vir¬ 
tually independent of the given domain and the number of nodes in the 
discretization of its boundary. In all numerical experiments we performed, 
we found that very few steps of (full) GMRES reduce the relative residual 
norm to 10 -12 or smaller. No preconditioning was required. Several exam¬ 
ples are given in Section [6] below. We believe that the very fast convergence 
of GMRES is due to the strong clustering of the eigenvalues of I — B around 
1 with only a few “outliers” and none of these “outliers” being close to zero. 
A more detailed analysis of this situation is a subject of further work. 
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In the numerical examples shown in Section [ 6 ] the MATLAB function 
fbie from [35] is used in order to obtain approximations to the unique 
solution fij of the integral equation (fTTl) and the function hj in (fT 2 l) . re¬ 
spectively. Within fbie we apply the MATLAB function gmres with the 
tolerance 10 ~ 14 for the relative residual norm. The matrix-vector product 
with I — B is computed using the MATLAB function zfmm 2 dpart from the 
MATLAB toolbox FMMLIB2D developed by Greengard and Gimbutas [TTI . 
We thus obtain approximations to the values of the functions pj and hj for 
j = !,...,£, at the points for i = 1 , 2, ... ,£n. Then the values of the 
constants hj are approximated by 

^ kn 

hj = ~ ^2 hh)- i = !,2, ...,£, k = 1,2, 

i=l-\-(k—l)n 

Since the function fbie requires 0(£n log n) operations, the computational 
cost for solving the l integral equations (fTTT) and computing the i functions 
h, in (1121) is 0(£ 2 nlogn) operations. For more details, we refer the reader 

to G3E51I3B!. 

Next, the values of the parameters mi ,..., mg, log t are computed by 
solving the linear algebraic system (1161) . Since £ usually is not large, this 
(£ + 1 ) x (£ + 1 ) system can be solved directly, e.g., using backslash in 
MATLAB. Finally, the values of the functions 7 and p at the points t t for 
i = 1 , 2, ..., £n can be computed from © and (fT3l) . 

5.2 Computation of the a,j and the boundary values of $ 

In the preceding section, we have computed the parameters m±,..., me, log t 
and the values of the functions /r and hj at the points t, for i = 1 , 2 ,..., in. 
In this section, we shall compute the values of a \,..., an and the values of the 
function <&(rj(t)) at the points tj for i = 1 , 2 ,... ,£n, by solving a non-linear 
system of equations. 

5.2.1 The non-linear system 

Let 


Wi = pi = logr + 7 (ti) + i/i(fi), i = 1, 2 ,..., £n. (28) 

The pi are known from Section 15.11 We will compute w\ , W 2 , ■ ■ ■, w^ n and 
a±, 02 , • ■., d£. We have from (1211) the following £n non-linear algebraic equa¬ 


ls 


tions in the in + i unknowns w\, W 2 , ■ ■ ■, wg n , a\,... ,ag 


^2 m i log 

3 =1 


v(ti) - a 3 


= Pi, i = 1,2,..., in. 


(29a) 


By discretizing the integral in (|22l) . we also have the following i non-linear 
algebraic equations in the in + i unknowns w\, W 2 , ■ ■ ■, wg n , a \,..., ag, 


1 


in 


— X! log 


3 = 1 


v(tj) - ai 


rj(tj) + ati — a,i = 0, i = l,2,...,£. (29b) 


Let z be the vector of the in + i unknowns w±, W 2 , ■ ■ ■, wg n , a \,..., at, 

i.e., 

Z = [wi ... w in ai ... ag] T G C £n+e , 
and let the function F be dehned by 


F( z) = 


E-=i^log(^ ] ^-) - Pi 


Ej=i lo § 

1 \^£n 


in Q'j 

Visin') 

Wj—ai 

j 


Pin 




G C 


£n+£ 


a Ej=i log (^r^) 

Then the system of non-linear equations (1291) can be written as 

F{ z) = 0. (30) 

5.2.2 Solving the non-linear system (1301) 

We shall solve the non-linear system (1301) using Newton’s iterative method 
z k+1 = z k - F\z k ) _1 F(z fc ), k = 0,1,2,... , (31) 

where F’( z) is the Jacobian matrix of the function F and is given by 


F’( z) = 


D Ai 

A-2 —Ig 


G C 


£n+£,Pn+£ 
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where is the t x l identity matrix, and 


D = 


mj 

2 —‘j=1 uh-a. 


E l TTlj 

7 = 1 W2~ai 


A, = 


Aq = 


zu 


W£n dj - 


6 C en ’ in , 


r —mi 

—m 2 

—m£ -| 


wi — ai 

Wl —CL2 

Wl —CL£ 


—mi 

—m 2 

—m£ 


W2—CLI 

W2—CL2 

W2 —CL£ 

£ C 

—mi 

— m 2 

—m£ 


“ w £n fil 

W£ n —CL2 

w en -a e -i 


ri f)(ti) 

1 »7(*2) 

1 v(ten) 1 

n\ w\ —a\ 

ni W 2 —ai 

ni W£ 


i fi(ti) 

i vfa) 

1 V{hn) 

n\ w\ —CL 2 

ni W 2 —CL 2 

ni W£ 

n~ a 2 

1 >?(ii) 

1 >j(*2) 

1 v(te n ) 

Lni w\ —d£ 

ni W 2 ~cl£ 

ni W£ 

n—a-e - 1 


'tin,l 


€ C 


l,In 


Let us discuss the choice of the starting point z°. In our numerical 
examples, a good choice for the starting point of a 3 has been found to be 
the center of mass of the boundary curve Tj, scaled by some factor > 1. A 
good choice for the starting point for the boundary values has been found 
to be small circles around aj. In the following we assume that a suitable 
starting point z° for the Newton method is used, so that, in particular, the 
matrix F'(z k ) is invertible in each iteration step. 

For each iteration k in (13111 . it is required to solve the linear system 


F'( z fc )v = F( z k ) 


(32) 


for v € C en+e . By taking into account the block structure of F’( z), the 
system (I32p can be reduced to an l x i linear system, as we show next. 

The vectors v and F(z k ) can be partitioned as 


v = 


F(z k ) = 


where x, b € C in and y,c £ C^. Hence equation (13211 is equivalent to the 
linear system 


£>x + Aiy = b, 
-42X - y = c. 


(33) 
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( 34 ) 


Lemma 5.1. The diagonal matrix D = \dij\ € <C £n,£n satisfies 

i ..-UAA j _ 12 tn 

u n — TT/ x j 1 — - 1 -; ‘‘i ■ ■ ■ i rn, 

U(Wi ) 

where U(w ) = f([y_i(u; — Pj) mj > the fi\, @ 2 , ■ ■ ■, fit are the last t entries of 
z k , and mi,... ,me are the exponents of the lemniscatic domain C. 

Proof. For the function U, we have 


U'(w) 

U(w) 


^ logU{w) 


_d_ 

dw 


l t 

- a i ) = X] 

j= 1 i =1 


m 3 

W — Qj 


Substituting w = Wi shows 4MD . □ 

In view of Lemma 12.21 the lemma suggests that D is non-singular if z k 
is close to the solution of F(z) = 0. This has also been observed in all 
numerical experiments; see Section [HI 

If D is not singular, we can rewrite the first equation in (|33l) as 

x = D -1 (b - Aiy) (35) 


and insert this in the second equation in fl33j) to obtain 

(A 2 D- 1 A 1 +I e )y = A 2 D- 1 b-c. (36) 


Now, we get a solution of fl32j) by solving the £xl system (136)) and computing 
x by (1351) . instead of solving the (in + £) x (in + i) system (1321) directly. 
The i x l system (l36j) can be solved using a direct method such as the 
Gauss elimination method since l is usually small. If D is non-singular, the 
matrix A 2 D l A\ + is invertible if and only if F'(z k ) is invertible. Thus 
the system (136T) is then uniquely solvable. 

Next we show that it is possible to compute the vectors x in (|35l) and y 
in (13611 without forming the matrices A\, A 2 or D first, by taking into account 
the Cauchy structure of A± and A 2 . Indeed, with the Cauchy matrix 


C = 


1 

1 

1 -| 

wi—ai 

1 

W1—CL2 

wi-a£ 

W 2 —CL\ 

W 2 —C 12 

W 2 —d£ 

1 

1 

1 

Win—a 1 

W£ n -a 2 

Wg n d £- 


= r—-—1 e c £n ’ e 

LWr CLs J vs ^ 
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the matrices A\ and A 2 can be written as 


777i 



777 2 

1 WT 

V(h) 


, a 2 — —C 



771 


m i. 


V(tin)_ 


Then the matrix A 2 D l A\ in the linear system (I36j) can be written as 


A 2 D~ 1 A l = - 
n 

so that we can generate the entries of this matrix directly, without first 
forming A\ or A 2 . Similarly we can write the right-hand side of (1361) as 


E tn _ ms _ 

k=l (w k -a r )(w k -a 3 ) -yt ln J 

'.7 = 1 Wb.—a a 


<E C* 


a 2 d 


-1 


b — c = — 
ni 


E 


with b = [ 61 , 62 , • • • , 6 te ] r . 
from ([35]) . 


in b k _ i)(t k ) 

k =1 w k —a r y-' '"j 

^3 = 1 w k~ a j 


r=1,2,.../ 


C£C‘ 


Finally, for computing the vector x, we have 


x = D 1 (b- Ai y) 


1 





m kVk 

W r Q>k 


r=l,2,...,£n 


withy = [yi,y 2 , ■ ■ ■, ye] T ■ 

The pseudo-code in Algorithm Q] summarizes our method described in 
Sections 15.1115.21 

We have used this method in the numerical experiments shown in Sec¬ 
tion [ 6 j 


5.3 Computation of the interior values of $ 


The method described above yields boundary values of the function <3?, 
namely the values <h( 7 /(t)) at the points tj for i = 1,2 ,...,£n. The val¬ 
ues of at interior points z € 1C can be computed by Cauchy’s integral 
formula applied to the function <l?(z) — z, which is analytic throughout /C 
and vanishes at 00 , 




77 (f) - z 


dt. 


A fast and accurate method to compute the Cauchy integral formula has 
been given in |35[ (see also 02 [18], [36]). The method is based on using the 
MATLAB function zfmm2dpart in [T3]. To compute the Cauchy integral 
formula at p interior points, the method requires 0 (p + £n ) operations. 
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Algorithm 1 Pseudo-code of the method 

Input: discretization t of J as in (1271) . parametrization 77 (f), 77 (f) of the 
boundary of /C, and auxiliary points a.j in the interior of the boundary 
curves Fj (j = 1 , 

Output: boundary values 3>(??(t)) and parameters 01,..., an, mi,.. . , m^, r 
of the lemniscatic domain. 

For j = 1,... ,£\ 

Compute 7 j(t) = — log ( 77 (f) — otj | from ©. 

Compute fij( t) and hj(t) from the boundary integral equa¬ 
tion CD and (I12p with the MATLAB function fbie. 

End 

Compute mi, ..., log r by solving the linear algebraic system (fT 6 l) . 
Compute the pi from (|28j) . where y(t) and /r(t) are given by (]9|) 
and CD- 

Compute the boundary values of <J> and ai, ..., ai by solving the non¬ 
linear system (1301) with Newton’s method. 


Remark 5.2. The Cauchy integral formula can also be used to compute the 
values of the inverse mapping for interior points w G C fl9 : p. 380]. 
However, it requires the boundary values of both the function and its 
derivative d*', i.e., &(r](t)) and &'(r](t)) at the points ti for i = 1,2,... ,ln. 
For computing the boundary values of <h / , we can use the boundary integral 
equation with the adjoint Neumann kernel as in m- Alternatively, we can 
compute Q'(r](t)) numerically from <&(r](t)). 

6 Numerical examples 

In this section we present numerical examples for five domains that illustrate 
our method. 

Example 6.1. We consider the unbounded domain /C exterior to the two 
circles 

771,2 (f) = ±1 + re” 1 *, 0 < t < 27r, 0 < r < 1, 

for different values of the radius r. The corresponding conformal map <I> and 
lemniscatic domain 

C = {wG C:\w- ai| mi \w - a 2 | m2 > r} 
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(both depending on the value of r ) have been derived analytically in |[42l Sec¬ 
tion 4], We therefore can compare our numerically computed parameters 
with the exact values aq, 02 , m±, m 2 , r defining C. Figure [2] shows the do¬ 
mains 1C for r = 0.5, r = 0.7, and r = 0.9 and the corresponding lemniscatic 
domains C. 

We denote the numerically computed parameters of the lemniscatic do¬ 
main by ai )n , a 2) n, m\, n ,m 2) n, T n , where n is the number of nodes in the 
discretization of each boundary component. The (absolute) errors 

Ba,n = max{|ai Ol,n|> |®2 
E m , n = max{|mi - mi, n |, \m 2 - m 2 , n |}, 

Er,n — \t T n \ 

are shown in Figure [3j We observe that all errors are quite small already 
for a small number of nodes. In fact, with n = 2 6 = 64 nodes the errors in 
this example are close to the machine precision level of 10~ 16 . Increasing 
the number of nodes beyond this point leads to some irregularities in the 
observed convergence behavior, which most likely is due to the some slight 
differences in the accuracy of the computed solution of the respective linear 
algebraic systems. Since all errors remain on the order of 10 -14 or smaller, 
we did not further investigate this phenomenon. 

Our method requires solving one linear algebraic system with the matrix 
I — B of size in x in with a different right hand side for each boundary 
component. In this example we have i = 2, and hence there are two linear 
algebraic systems to be solved for each value of r and n. As described in 
Section 15.11 we use the (full and unpreconditioned) GMRES method for this 
task. In Figure [4] we plot all relative residual norms of the GMRES method 
we obtained for the two linear algebraic systems, r = 0.5, 0.7,0.9, and n = 
2 6 ,2 7 ,..., 2 10 . Thus, Figure Q] shows the GMRES convergence for 30 linear 
algebraic systems. We observe that the number of GMRES iteration steps 
required to attain a relative residual norm on the order of 10 -14 is very 
small and almost independent of parameters in the linear algebraic systems 
(namely the right hand side, r, and n). We have indicated reasons for this 
observation in Section ECU but, as mentioned there, a detailed analysis is 
the subject of future work. 

The 2-norm condition numbers of the matrices D and A 2 D~ l A\ + / (_ in 
the Newton iteration are shown in Figure 0 All matrices stay quite well 
conditioned throughout the iteration. The same observation can be made 
in the following examples as well. Figure [6] shows the norms ||z fc+1 — z fc || 00 
in the Newton iteration. 
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Figure 2: Original domains for Example [6J] (left) and corresponding lemnis- 
catic domains (right) obtained with n = 256 and for r = 0.5 (top), r = 0.7 
(middle), and r = 0.9 (bottom). 
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Figure 3: Errors E a>n , E m _ n , and E r ^ n obtained with n nodes and for r = 0.5 
(top left), r = 0.7 (top right), and r = 0.9 (bottom) in Examp le i 6. II Missing 
dots indicate that the error is zero. 



Figure 4: Relative residual norms of the GMRES method for different values 
of r and n in Example 16. II 
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Condition number 



Number of iteration: k 

Figure 5: 2-norm condition numbers of D and A 2 D~ l A\ + Ig obtained with 
n = 256 and for r = 0.5 (top left), r = 0.7 (top right), and r = 0.9 (bottom) 
in Example 16.11 



Figure 6: ||z fc+1 — in the Newton iteration obtained with n = 256 and 
for r = 0.5, r = 0.7, and r = 0.9 in Example 16.11 
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Figure 7: Original domain for Example 16.21 (top) and corresponding lemnis- 
catic domain (bottom) obtained with n = 256. 


Example 6.2. We consider the unbounded domain /C exterior to seven 
nonconvex and complicated but smooth curves as shown in Figure [71 These 
curves are parametrized (from left to right) by 

Vj(t) = rj(t) e _lt , 0<f<2vr, j = 1,2,..., 7, 


where 


r\(t ) = 1.25 + 0.50sin(4f) + 0.30cos(t), 

7*2 (t) = 1.25 + 0.40 sin(2t) + 0.20 cos(3i), 

r 3 (t) = 0.75 + 0.25e cos(i) cos 2 (3f) + 0.50e sin(t) sin 2 (2f), 

r 4 (t) = e cos W cos 2 (2t) + e sin W sin 2 (2f), 

r 5 (t) = 0.75 + 0.25e cos(t) cos 2 (2t) + 0.50e sin(t) sin 2 (3t), 

re(t) = 1.25 + 0.40sin(4t) + 0.20cos(3t), 

ri{t) = 1.25 + 0.50sin(3f) + 0.30cos(t). 

The computed lemniscatic domain obtained with n = 256 is shown on the 
bottom of Figure [TJ The GMRES method for the seven linear algebraic 
systems required 36 iteration steps to attain a residual norm smaller than 
10 -14 . Figure [8] shows the 2-norm condition numbers of the matrices D and 
A 2 D~ l A\ + 1£ as well as the norms ||z^' +1 — z fc || 00 in the Newton iteration. 

Example 6.3. In this example we demonstrate that our method also works 
for domains with high connectivity. We consider the unbounded domain /C 
exterior to 64 circles as shown on the left of Figure [9j The domain /C is 
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Figure 8: 2-norm condition numbers of D and A 2 D 1 A\ + In (left) and 
norms ||z fc+1 — z^Hqo (right) in Example 16.21 


symmetric with respect to both the real and the imaginary axis. Since the 
computed conformal map is normalized as in ([3]) the lemniscatic domain 
C has the same symmetry properties; see [121 Lemma 2.2]. The computed 
lemniscatic domain obtained with n = 256 is shown on the right of that 
figure. The GMRES method for the 64 linear algebraic systems required 
between 14 and 18 iteration steps to attain a residual norm smaller than 2 ■ 
10 -14 . Figure [10] shows the 2-norm condition numbers of D and A?D~ 1 A\ + 
In as well as the norms ||z^ +1 — z fc ||oo in the Newton iteration. 

Example 6.4. As indicated in Section 15.11 our method can also be used 
when the boundary components of /C are only piecewise smooth Jordan 
curves. As an example we consider the unbounded domain /C exterior to 
four squares as shown on the left of Figure [TTJ The computed lemniscatic 
domain obtained with n = 1024 is shown on the right of that figure. The 
GMRES method for the four linear algebraic systems required 26 iteration 
steps to attain a residual norm smaller than 10~ 14 . Figure [T2l shows the 
2-norm condition numbers of D and A 2 D~ 1 A\ + In as well as the norms 
||z^ -1-1 — z fc || 00 in the Newton iteration. 

Example 6.5. In this example we consider the unbounded domain 1C exte¬ 
rior to three non-convex sets as shown on the left of Figure [T] The three sets 
are of the form introduced in [22]. The boundary curves are analytic and an 
analytic parameterization is known as well. The corresponding lemniscatic 
domain obtained with n = 1024 nodes per boundary component is shown 
on the right of Figure [U The GMRES method for the three linear algebraic 
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Figure 9: Original domain for Example 16.31 (left) and the corresponding 
lemniscatic domain (right) obtained with n = 256. 



Figure 10: 2-norm condition numbers of D and A 2 D l A\ +Ii (left) and the 
norms — z^H^ (right) in Example 16.31 
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Figure 11: Original domain for Example 16.41 Cleft) and corresponding lem- 
niscatic domain (right) obtained with n = 1024. 




Figure 12: 2-norm condition numbers of D and A 2 D l A\ +1% (left) and the 
norms (right) in Example 16.41 
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Figure 13: 2-norm condition numbers of D and A 2 D l A\ +Ig (left) and the 
norms ||z fc+1 — z fc || 00 (right) in Example 16.51 


systems required between 32 and 34 iteration steps to attain a residual norm 
smaller than 10~ 14 . Figure fT3l shows the 2-norm condition numbers of D and 
A 2 D~ x A\ + Ip as well as the norms Hz^ -1-1 — z*'||oo in the Newton iteration. 


7 Concluding remarks 

In this article we derived a method that numerically computes the conformal 
map from a given domain onto a lemniscatic domain. The method relies on 
solving a boundary integral equation with the Neumann kernel. It takes 
as input a parameterization of the boundary of the original domain and 
yields the parameters of the lemniscatic domain and the boundary values 
of the conformal map. Using the numerically computed conformal map : 
/C —y T it is in particular possible to compute the Faber-Walsh polynomials 
associated with the compact set C \ /C. The first numerical examples for 
such a computation are given in the paper [43l . 

The transfinite diameter (or logarithmic capacity) of the set C\/C is one of 
the parameters in its corresponding lemniscatic domain £; see Theorem 12.11 
This parameter is computed in the first step of the algorithm proposed in 
this paper; see Section 15.11 Since computing the transfinite diameter of 
compact sets is an interesting problem in its own right, we have derived a 
numerical method for this task, which is based on the method of Section [5.11 
in our paper m- 

Let us point out a few open questions. As discussed in Section 15.11 
and demonstrated numerically in Section [6l the GMRES method converges 
very fast when solving the discretized boundary integral equation with the 
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Neumann kernel. A rigorous analysis of this effect is subject to further work. 
Further, it would be interesting to analyze how the accuracy of the solution 
of the linear algebraic systems (solved with GMRES) affects the accuracy 
of the computed conformal map and of the parameters of the lemniscatic 
domain. Finally, a “black box” starting point for the Newton iteration for 
solving the non-linear system (1301) would be of interest. 
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